Retirado do artigo Miller et al. (2019). Distance sampling in R. Journal of Statistical Sofware 89(1)
# instalar pacotes necessários
#install.packages("Distance")
# instalar pacotes adicionais
#install.packages("mrds")
#install.packages("dsm")
#install.packages("mads")
#install.packages("dsims")
# carregar pacotes
library(Distance)
Carregando pacotes exigidos: mrds
This is mrds 2.2.8
Built: R 4.3.0; ; 2023-04-28 17:39:52 UTC; unix
Attaching package: ‘Distance’
The following object is masked from ‘package:mrds’:
create.bins
library(dplyr)
library(DT)
library(flextable)
Attaching package: ‘flextable’
The following object is masked from ‘package:purrr’:
compose
The following objects are masked from ‘package:plotly’:
highlight, style
The following objects are masked from ‘package:ggpubr’:
border, font, rotate
library(ggplot2)
library(lubridate)
library(plotly)
library(readr)
library(readxl)
library(stringr)
library(tibble)
library(tidyr)
# carregar as funções da pasta R
# carregar função script_carregar_funções_pasta_r.R
source(
paste0(
here::here(),
"/R/minhas_funcoes.R"
)
)
# carregar dados
cutia_tap_arap <- transformar_para_distanceR_covariaveis() |>
filter(
uc_name == "Resex Tapajos-Arapiuns",
sp_name == "Dasyprocta croconota"
) |>
drop_na(distance)
# readr::write_excel_csv(
# cutia_tap_arap,
# paste0(
# here::here(),
# "/data/cutia_tap_arap.csv"
# ),
# )
cutia_tap_arap |>
DT::datatable(filter = "top")
cutia_tap_arap |>
count(Region.Label)
cutia_tap_arap |>
count(Sample.Label)
Variáveis necessárias para o data.frame:
Region.Label: vetor fator com o estrato contendo o
transecto (pode ser uma estratificação pré-amostragem - UCs - ou
pós-amostragem - ex. região, estado, bioma)
Area: vetor numérico contendo a área do
estrato;
Sample.Label: vetor númerico contendo a identidade
(ID) do transecto
object: nome adicional, ver seção 6;
detected: nome adicional, ver seção 6;
Effort: vetor númerico contendo o esforço do
transecto (para linhas seu comprimento, para pontos o número de vezes
que o ponto foi visitado)
size: vetor numérico copntendo o tamanho do grupo
observado;
distance: vetor numérico de distâncias
observadas;
Month:
OBs:
Sp:
mas:
HAS:
Study.Area:
Transectos que foram amostrados, mas que não tiveram observações (n =
0) devem ser incluídos no conjunto de dados com NA nas
observações de distância e qualquer outra covariael para a qual não se
tenha observação.
# cutia_tap_arap |>
# complete(Region.Label, Sample.Label, sp_name) |>
# datatable(filter = list(position = "top"))
Jogar a imputacao de NAs pra dentro da funcao carregar
dados completos.
# desenha o grafico com a distribuicao de distancias perpendiculares
cutia_tap_arap |>
plotar_distribuicao_distancia_interativo()
Warning: Continuous y aesthetic
ℹ did you forget `aes(group = ...)`?
Ajustando um modelo ao dados das cutias Dasyprocta
croconota, configurando uma distância limite de 20m e usando
Half-normal como key function usando o argumento
key, sem termo de ajuste.
# ajustando a função de detecção para uma distancia de truncamento de 20 metros
# Key function - Half-normal
cutia_tap_arap_hn <- cutia_tap_arap |>
ds(
truncation = 20,
key = "hn",
adjustment = NULL
)
Fitting half-normal key function
AIC= 7227.642
Warning: Some observations not included in the analysis
Ajustando um modelo ao dados das cutias Dasyprocta
croconota, configurando uma distância limite de 20m e usando
Half-normal como key function usando o argumento
key, com termo de ajuste do tipo
Hermite-polynomial.
# ajustando a função de detecção para uma distancia de truncamento de 20 metros
# Key function - Half-normal
cutia_tap_arap_hn_herm <- cutia_tap_arap |>
ds(
truncation = 20,
key = "hn",
adjustment = "herm"
)
Starting AIC adjustment term selection.
Fitting half-normal key function
AIC= 7227.642
Fitting half-normal key function with Hermite(4) adjustments
AIC= 7229.142
Half-normal key function selected.
Warning: Some observations not included in the analysis
Ajustando um modelo ao dados das cutias Dasyprocta
croconota, configurando uma distância limite de 20m e usando
Half-normal como key function usando o argumento
key, com termo de ajuste do tipo Cosine.
# ajustando a função de detecção para uma distancia de truncamento de 20 metros
# Key function - Half-normal
cutia_tap_arap_hn_cos <- cutia_tap_arap |>
ds(
truncation = 20,
key = "hn"
)
Starting AIC adjustment term selection.
Fitting half-normal key function
AIC= 7227.642
Fitting half-normal key function with cosine(2) adjustments
AIC= 7207.469
Fitting half-normal key function with cosine(2,3) adjustments
AIC= 7201.248
Fitting half-normal key function with cosine(2,3,4) adjustments
Warning: Detection function is not strictly monotonic!AIC= 7152.177
Fitting half-normal key function with cosine(2,3,4,5) adjustments
Warning: Detection function is not strictly monotonic!AIC= 7146.681
Fitting half-normal key function with cosine(2,3,4,5,6) adjustments
Warning: Detection function is not strictly monotonic!AIC= 7129.427
Warning: Detection function is not strictly monotonic!Warning: Some observations not included in the analysis
Ajustando um modelo ao dados da cutia Dasyprocta croconota,
configurando uma distância limite de 20m e usando Hazard rate
como key function usando o argumento key.
# Key function - Hazard-rate
cutia_tap_arap_hr <- cutia_tap_arap |>
ds(
truncation = 20,
key = "hr",
adjustment = NULL
)
Fitting hazard-rate key function
Warning: Estimated hazard-rate scale parameter close to 0 (on log scale). Possible problem in data (e.g., spike near zero distance).AIC= 6907.541
Warning: Estimated hazard-rate scale parameter close to 0 (on log scale). Possible problem in data (e.g., spike near zero distance).Warning: Some observations not included in the analysis
Ajustando um modelo ao dados da cutia Dasyprocta croconota,
configurando uma distância limite de 20m e usando Hazard rate
como key function usando o argumento key e termo
de ajuste do tipo polinomial simples.
# Key function - Hazard-rate
cutia_tap_arap_hr_poly <- cutia_tap_arap |>
ds(
truncation = 20,
key = "hr",
adjustment = "poly"
)
Starting AIC adjustment term selection.
Fitting hazard-rate key function
Warning: Estimated hazard-rate scale parameter close to 0 (on log scale). Possible problem in data (e.g., spike near zero distance).AIC= 6907.541
Fitting hazard-rate key function with simple polynomial(4) adjustments
Warning: Estimated hazard-rate scale parameter close to 0 (on log scale). Possible problem in data (e.g., spike near zero distance).Warning: Model fitting did not converge. Try different initial values or different model
Model failed to converge.
Hazard-rate key function selected.
Warning: Estimated hazard-rate scale parameter close to 0 (on log scale). Possible problem in data (e.g., spike near zero distance).Warning: Some observations not included in the analysis
Ajustando um modelo ao dados da cutia Dasyprocta croconota,
configurando uma distância limite de 20m e usando Hazard rate
como key function usando o argumento key e termo
de ajuste do tipo cosseno.
# Key function - Hazard-rate
cutia_tap_arap_hr_cos <- cutia_tap_arap |>
ds(
truncation = 20,
key = "hr"
)
Starting AIC adjustment term selection.
Fitting hazard-rate key function
Warning: Estimated hazard-rate scale parameter close to 0 (on log scale). Possible problem in data (e.g., spike near zero distance).AIC= 6907.541
Fitting hazard-rate key function with cosine(2) adjustments
Warning: Estimated hazard-rate scale parameter close to 0 (on log scale). Possible problem in data (e.g., spike near zero distance).Warning: Model fitting did not converge. Try different initial values or different model
Model failed to converge.
Hazard-rate key function selected.
Warning: Estimated hazard-rate scale parameter close to 0 (on log scale). Possible problem in data (e.g., spike near zero distance).Warning: Some observations not included in the analysis
plot(cutia_tap_arap_hn, breaks = seq(0, 20, 2.5))
plot(cutia_tap_arap_hn_herm, breaks = seq(0, 20, 2.5))
plot(cutia_tap_arap_hn_cos, breaks = seq(0, 20, 2.5))
plot(cutia_tap_arap_hr, breaks = seq(0, 20, 2.5))
plot(cutia_tap_arap_hr_poly, breaks = seq(0, 20, 2.5))
plot(cutia_tap_arap_hr_cos, breaks = seq(0, 20, 2.5))
Ajustando um modelo ao dados das cutias, configurando uma distância
limite de 20m e usando Uniform como key function a
séria de ajuste com o argumento adjustment ee especificando
a ordem com o argumento order.
cutia_unifcos <- cutia_tap_arap |>
ds(truncation = 20,
key = "unif",
adjustment = "cos",
order = c(1, 2))
Ajuste Hermite pollynomial usa od código "herm"
e polinomial simples "poly".
Podemos incluir covariáveis utilizando o argumento
formula = ~ .... Abaixo, está especificado um modelo
“Hazard-rate” para os dados de cutia q ue inclui o tempo de senso como
covariável e uma distância limite de 20m.
cutia_hr_time <- cutia_tap_arap_15 |>
ds(truncation = 20,
key = "hr",
formula = ~ cense_time)
Adicionando uma segunda covariável: tamanho do grupo.
cutia_hr_time_size <- ds(data = cutia_tap_arap_15,
truncation = 20,
transect = "line",
key = "hr",
formula = ~ cense_time + size)
plot(cutia_hr_time)
plot(cutia_hr_time_size)
Podemos usar a função summary para obter informações
importantes sobre o modelo.
summary(cutia_hn)
O resultado inclui detalhes sobre o dado e a especificação do modelo, assim como dos coeficientes (\(\beta_{j}\)) e sua inceteza, a média do valor de detectabilidade e sua incerteza e uma estimativa da abundância na área coberta pela amostragem (sem levar em consideração o tamanho dos agrupamentos, ou bandos).
Para visualizar quão bem a função de detecção se ajusta aos dados quanto temos as distâncias exatas podemos usar um plot de quantis empíricos x teóricos (Q-Q plot). Ele compara a função de distribuição cumulativa (CDF) dos valores ajustados da função detecção a distribuição empírica dos dados (EDF).
Também podemos usar o teste de Cramér-von Mises para testar se os pontos da EDF e da CDF tem origem na mesma distribuição. O teste usa a soma de todas as distâncias entre um ponto e a linha y = x para formar a estatística a ser testada. Um resultado significativo fornece evidência contra a hiipótese nula, sugerindo que o modelo não se ajusta bem aos dados.
# ajustando um modelo Half-normal
cutia_hn <- ds(data = cutia_tap_arap_15,
truncation = 20,
transect = "line",
key = "hn",
adjustment = NULL)
# conduzindo o teste dfe bondadede ajuste de Cramer-von Mises
gof_ds(cutia_hn)
gof_ds(cutia_hr_time)
O resutlado do teste aponta que o modelo Half-normal deve ser descartado.
Testes de bondade de ajuste de chi-quadrado são gerados usando a
função gof_ds quando as distâncias forneceidas estão
categorizadas.
Uma vez que temos um conjunto de modelos plausíveis, podemos utilizar
o cirtériode informaçãode Akaike (AIC) para selecionar entre os modelos
o que melhor se ajusta aos dados utilizando a função
summarize_ds_models.
# gerando uma tabela de seleção de modelos usando AIC
summarize_ds_models(cutia_hn, cutia_hr_time, cutia_hr_time_size)
O melhor modelo é o Hazard-rate com tempo de senso e tamanho do grupo como covariáveis.
Para obter a abundância na região de estudo, primeiro calculamos a abundância na área amostrada para obter \(N_c\) e em seguida escalonamos esse valor para toda a área de estudo multiplicando \(N_c\) pela razão entre a área amostrada e a área da região. Para estimar a abundância na área amostrada, utilizamos as estimativas de probabilidade de detecção no estimador de Horvitz-Thompson.
Quando fornecemos os dados no formato correto (“flatfile”)
ds irá automaticamente calcular as estimativas de
abundância baseado nas informações de amostragem presenta nos dados.
summary(cutia_hn)
Summary statistics: fornece as áreas, aŕea de amostragem, esforço, número de observações, número de transectos, taxa de encontro, seus erros padrões e coeficientes de variação para cada estrato;
Abundance: fornece estimativas, erros padrões, coeficientesde variação, intervalos de confiança inferior e superior, graus de liberdade para a estimativa de abundância de cada estrato;
Densidade: lista as mesmas estatísticas de Abundance, só que para densidade.
contar_n_repeticoes_trilha() - conta o número de vezes
que cada trilha foi visitada